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Abstract 

Of several iterative and direct equation solvers evaluated previously for computa- 
tions in aeroacoustics, the most promising was the NASA-developed General-Purpose 
Solver (winner of NASA’s 1999 software of the year award). This paper presents de- 
tailed, single-processor statistics of the performance of this solver, which has been 
tailored and optimized for large-scale aeroacoustic computations. The statistics, com- 
piled using an SGI ORIGIN 2000 computer with 12 Gb available memory (RAM) and 
eight available processors, are the central processing unit time, RAM requirements, 
and solution error. The equation solver is capable of solving 10 thousand complex 
unknowns in as little as 0.01 sec using 0.02 Gb RAM, and 8.4 million complex un- 
knowns in slightly less than 3 hours using all 12 Gb. This latter solution is the largest 
aeroacoustics problem solved to date with this technique. The study was unable to 
detect any noticeable error in the solution, since noise levels predicted from these solu- 
tion vectors are in excellent agreement with the noise levels computed from the exact 
solution. The equation solver provides a means for obtaining numerical solutions to 
aeroacoustics problems in three dimensions. 
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1 Introduction 

Fan noise is defined as unwanted acoustic energy that is generated at the fan face or guide 
vanes of a turbofan engine. Fan noise accounts for a significant portion of community noise 
radiated from conventional and high bypass ratio engines. Noise reduction research today 
focuses on reducing fan noise levels radiated from future aircraft by a factor of 2 relative 
to 1990 levels. Installation of acoustic treatment (i.e., liners) into the nacelles of aircraft 
engines remains one of the most effective means for achieving these noise reduction goals [1], 
To this end, an accurate knowledge of liner impedance is critical in optimizing the treatment 
for maximum noise suppression. To date, much of the design effort has concentrated on 
expensive and time consuming experimental testing. In addition, experimental tests have not 
account for variable surface impedance that results naturally from imperfections in the liner 
manufacturing process. An accurate numerical model could predict the lining impedance 
in a less costly and more time-efficient manner, and at the same time account for surface 
impedance variability. 

In a recent paper [2] a numerical method for extracting the impedance of an acoustic 
material was developed and validated for plane wave sources. In this approach, the time- 
dependent acoustic equations are Fourier transformed into a single differential equation in 
frequency space. Source and exit boundary conditions for the numerical model are obtained 
from measurements in a flow duct that provides grazing-flow and grazing-incidence sound 
over the test liner. The frequency domain differential equation is then coupled to the mea- 
sured boundary conditions and the solution is approximated by a standard finite element 
method. The finite element method leads to a large, sparse, indefinite linear system of 
complex equations. The acoustic impedance of the test liner is then educed by a series of 
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successive corrections to an original estimate for the unknown impedance, the process being 
carried out repetitively until the corrected impedance reproduces the measured upper wall 
pressure. Each successive correction to the impedance is determined via an optimizer that 
uses the solution to the system of complex equations to construct the penalty function [2]. 
As this method has progressed to include higher frequencies, nonplanar sound sources and 
flow [3], the absence of an efficient (in time and memory) equation solver for this large system 
of complex acoustic equations has emerged as a major impediment to further development 
of the method. A time-and memory-efficient solver for aeroacoustics would allow numeri- 
cal methods to be extended to high frequency sources and three-dimensional aeroacoustics 
computations. 

Large systems of complex equations may be solved by iterative [4] or direct [5] methods. 
Recently, the performance of several iterative and direct equation solvers were evaluated to 
establish their suitability for computations in aeroacoustics [6]. Based on that study, the 
commercial version of the NASA-developed General-Purpose Solver (GPS) [7] emerged as 
the most promising solver. However, the conclusion derived in [6] was based solely on a 
study of the central processing unit time required by the solver. Among the other metrics 
that need to be considered in evaluation of the solver are memory requirements (RAM) 
and solution error. Additionally, for realistic aeroacoustics computations, the impedance 
spectrum of a test sample would be required for frequencies up to 25,000 Hz and 0.5 Mach 
number. A “back of the envelope calculation” at this frequency and Mach number shows 
that approximately 8 million equations need to be solved to resolve all propagating modes 
contained within the computational volume. Unfortunately, due to memory constraints the 
results presented in [6] were only small-scale calculations (i.e., less than 80,000 equations), 
and no attempt was made to obtain optimal statistics. 

The purpose of this paper is to present detailed statistics of the performance of a vastly 
improved version of the GPS solver [8] that has been specifically tailored and optimized for 
nacelle aeroacoustics computations. The three metrics used in this evaluation are CPU time, 
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RAM requirements, and solution error. These statistics are evaluated for unknowns ranging 
from 10 thousand to 8.4 million. 


2 Physical Problem Description 


Figure 1 shows a schematic of the two-dimensional duct used in this study. It should be noted, 
as suggested by Figure 1, that the math model discussed here is limited to a two-dimensional 
description that approximates a three-dimensional flow impedanc# tube as discussed in [3] . 
The axial and transverse directions are denoted by .r and y, respectively. The duct is L units 
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Fig. 1. Two-dimensional duct and coordinate system. 


long, with the source and exit planes located at x = 0 and x = L, respectively. Required 
inputs at the source and exit planes are the source pressure p s {y) and the normalized exit 
impedance Cexit(y), respectively. Throughout this work all impedances are normalized with 
respect to the characteristic impedance, p 0 Qn °f the air in the duct, and the upper wall is 
rigid. Here p 0 and c 0 are the mean density and sound speed of air respectively, in the duct. 
Note that there is a mean flow in the axial direction that flows subsonically from left to right 
with uniform speed, u 0 . Further, there are m points, located at x = .iq, x- 2 , x 3 , . . . , x m along 
the upper wall, at which the acoustic pressures, p(xi,H) are assumed known. It should be 
noted that flow impedance tube apparatuses such as that discussed in [3] are routinely used 
to obtain measurement of the source pressure p s (y ), exit impedance Cexit(y), and the upper 
wall acoustic pressures, p(xi,H), for impedance;: eductions. 

The sound-absorbing material constitutes the part of the bottom wall of the duct between 
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L\ < x < L 2 . Outside this region, the lower wall is rigid. The acoustic material is assumed 
to be locally reacting (i.e., acoustic waves propagate through it only normal to the surface). 
The sound-absorbing material has an unknown normalized impedance ((x), as shown. The 
physical problem is to determine, from the measured boundary data, the impedance of the 
acoustic material as a function of the mean flow speed u 0 . 

3 Governing Equations and Boundary Conditions 

The equations that describe the propagation of acoustic pressure disturbances through a 
duct containing a flowing fluid as depicted in Figure 1 are derived from the Navier-Stokes 
and energy equation, neglecting viscous and heat-conducting effects. The justification for 
the neglect of viscosity and heat conduction is that the passage of sound waves through a 
moving fluid is an isentropic process. Note that the equations that are the subject of this 
investigation result from the additional assumptions that 

1. nonlinear acoustic effects can be neglected 

2. the acoustic disturbance has reached a steady-state 

3. the steady-flow is in uniform motion 

Thus, the mathematical problem is to find the solution to the convected wave equation [9] 

, = v VT; k = tL. M 0 = — (2) 

Co Co 

where p(x, y) is the complex acoustic pressure in the duct and / is the sound source frequency 
in hertz. The source, exit, and upper wall boundary conditions are [3] 

p{Q,y)=Ps{y) (3) 

9p(L,y) = ikp(L,y ) 

9x [Mo + Cexit(y)] 
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^ = 0 (5, 

while the correct form of the wall impedance boundary condition in the presence of the 
flowing fluid has been derived in [10] 
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Together, Equation (1) and the boundary conditions, Equations (3)-(6), form a boundary 
value problem that can be solved to determine uniquely the upper wall pressures p(xi,H) 
for a given impedance function. Conversely, if the upper wall pressures p(xi, H), are known 
(i.e. , measured), there is a unique wall impedance function £(a:) that will reproduce these 
wall pressures. This impedance function is the unknown impedance of the test liner. 

Unfortunately, exact solutions for the wall impedance function satisfying the boundary 
value problem do not exist for an arbitrary set of boundary data, p s (y), Cexit (y), and p(x, H). 
Thus, the goal of impedance eduction techniques is to devise a numerical procedure for 
determining this unknown liner impedance in the presence of the flowing fluid. In a re- 
cently developed method [2, 3], the unknown acoustic impedance function £(a:) of the test 
liner is “educed” by a series of successive corrections to an original estimate for the un- 
known impedance, the process being carried out repetitively until the corrected impedance 
reproduces the known upper wall pressure. Each successive correction to the impedance is 
determined by repeated numerical solutions to the nacelle acoustics boundary value problem 
defined by Equation (1) and boundary condition Equations (3)-(6). 

The process of obtaining the numerical solution to the boundary value problem consists 
of two stages. In the first stage the continuous partial differential equation and boundary 
conditions are converted into a system of discrete algebraic equations. The second stage re- 
quires an equation solver to obtain the solution to the system of discrete algebraic equations. 
A number of methods are available for converting the continuous differential equations and 
boundary conditions into a system of discrete algebraic equations. The method used here is 
the finite element method. 


6 



4 Formulation of Discrete Equation System 


The continuous partial differential equation and boundary conditions are converted into 
a system of discrete algebraic equations using the finite element method. As shown in 
Figure 2, N and M points are used in the x and y directions, respectively. Note that the 
computational domain has been decomposed into a total of (IV — l)x(M — 1) rectangular 
elements, as shown in the figure. Linear and cubic Hermite polynomial basis functions are 
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Fig. 2. Finite element discretization of two-dimensional duct. 


used with and without flow, respectively. Implementations of the finite element method are 
given in detail elsewhere [2, 3] and are not discussed further in this paper. The finite element 
method leads to a discrete set of complex linear equations of the form 


[A(C(r))]{$} = {F} 


( 7 ) 


When the duct does not contain flow (i.e., M 0 = 0), the coefficient matrix is of the form 


[A(C(*))] = 


[Ai] [B 2 \ 

[B 2 ] t [A 2 ] [B 3 ] 

[B N - i] t [Ajv_i] [B n ] 
[B n ] t [Ajv] 


( 8 ) 
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where the superscript T denotes matrix transpose and 


[Bi\ = 
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where aj, 6 j,Cj,dj 5 an( i e j are complex coefficients. Thus, the zero flow matrix [A(£(a:))], 
has the following properties 

1. It is an (N M)x(N M) indefinite matrix (i.e., construction of [L\ and [U] requires row 
or column interchanges in [A]) 


2. It is complex and symmetric 

3. Each major block ([A/], [Bj]) is an MxM complex tridiagonal matrix and the diagonal 
blocks ([A/]) are symmetric 


4. It is banded with a bandwidth of (2 M + 3) 

5. It is very sparse (i.e., only four of the (M + 1) superdiagonals contain nonzero elements) 
When the duct does contain flow (i.e., M 0 ^ 0), the coefficient matrix is of the form 
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where and [cj\ are 4x4 complex matrices. When the nacelle contains flow, the 

structure of [Aj] and \Ci\ are identical to that of [-B/], but the complex coefficients in [A/] 
and [Ci\ are distinct from those of [Bj], Thus, the flow matrix [A(£(x))], has the following 
properties 

1. It is an (4N M)x(4N M) indefinite matrix 

2. It is complex and nonsymmetric 

3. Each major block is a 4Mx4M complex matrix that is block tridiagonal 

4. It is banded with a bandwidth of 8 (M + 2) — 1 

5. It is very sparse 

Much practical importance arises from the structure of [A(£(x))], as it is convenient for 
minimizing storage and maximizing computational efficiency. Approximately 98% of the 
computer resources required to educe the impedance are consumed in Ending the solution 
to the discrete system. Considerable savings in computer resources are possible if the most 
efficient solver is used. In [2] and [3], the solution to the discrete system is obtained by using 
a band solver. However, the band solver severely taxes RAM and CPU time by requiring 
storage and arithmetic operations on the inner null bands of [A(£(x))]. This paper focuses on 
the use of an efficient sparse solver to obtain the solution of the system of discrete equations. 

5 The Equation Solver 

In this research effort, VSS [7] (a commercial version of NASA’s GPS [8]) is exercised to 
obtain the solution to the aeroacoustic system defined by Equation (7). The NASA-developed 
GPS had its genesis in the need for large aerospace structures solutions in computational 
mechanics. GPS was subsequently extended to support matrices that are sparse or dense, 
positive definite or indefinite, real or complex. In addition, GPS has been extended to solve 
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nonsymmetric matrices that are often generated, for example, in aeroacoustics problems 
containing flow. 

The method of solution used by GPS and VSS is the commonly-used [T][U] decomposi- 
tion. A fairly general scheme exists for computing the lower triangular factorization matrix 
[L\ and the upper triangular factorization matrix [U]. One key innovation of the GPS soft- 
ware is that [L\ and [U] are computed cleverly, with due attention given to eliminating 
computations with zero elements while minimizing storage and CPU time. A second GPS 
innovation is a novel reordering method that retains the benefit of a multiple-minimum de- 
gree (MMD) reordering at a fraction of the MMD reordering time. This time reduction is 
accomplished by reordering a subset of the equations. The solver requires that only the 
nonzero coefficients in [A(£(:r))] be stored in memory. The nonzero coefficients are stored 
in row format and as a single vector to facilitate the solution procedure. The GPS (and 
VSS) equation solver exploit the matrix characteristics (real or complex, symmetric or non- 
symmetric, in-core or out-of-core) of the application and is designed to exploit the hardware 
features of current and future computers. Only a small fraction of the capability of the 
solver is used in this research effort (i.e., only the complex, symmetric, in-core capability 
was required). The equation solver used in this paper includes several recent innovations 
that are discussed in detail elsewhere [8]. The statistics computed in the following section 
are optimal in the sense that several of the solver software parameters, such as the “loop 
level” [7], have been optimized through numerical experimentation. 

6 Results 

An in-house computer code that constructs the coefficient matrix [A(£(:r))] in the required 
NASA solver format has been linked to a commercial version of the GPS solver (VSS) in 
order to provide statistics for the aeroacoustics computations presented in this paper. The 
statistics for the solution of Equation (7) have been computed for both rigid-wall and soft-wall 
ducts, but in the absence of flow. Results were computed on an in-house SGI ORIGIN 2000 


10 



computer that contained 12 Gb of available RAM and eight processors. All computations 
were performed on a single processor using double-precision (i.e., 64 bit) arithmetic with 
M = N. 

Results are presented for a square duct 1 m in width (L = H =1 m) with the leading and 
trailing edges of the liner located at the source and exit planes, respectively = 0, L 2 = L). 
Because the results computed here do not contain flow effects (i.e., M 0 = 0), only the 
symmetric version of the solver is used. All calculations presented in this paper are performed 
at standard atmospheric conditions. The source is a plane wave source ( p s (y ) = 1) that 
oscillates at a frequency of 1000 Hz and the exit impedance is set to unity (Cexit = 1)- 
Because results in this paper have been purposely restricted to this range of parameters, 
exact solutions are available to check the solver error for the rigid- wall duct. 

Figure 3 shows a plot of the CPU time and RAM required to solve for up to 8.4 million 
complex unknowns in a rigid-wall duct. Note that these statistics are plotted on a dual axis 
system with the CPU time referenced to the y 1 axis and the memory referenced to the y 2 axis. 
CPU times shown in the figure correspond to those required to obtain the solution vector 
(i.e., to forward and backward solve) after the coefficient matrix is constructed. Generally, 
the CPU time required for construction of the coefficient matrix or to obtain the “backward 
solve” is less than 2% of that required to obtain the solution vector. CPU times range from 
0.01 sec for 10 thousand complex unknowns to a maximum of nearly 3 hours for 8.4 million 
complex unknowns. The RAM ranges from 0.02 Gb for 10 thousand complex unknown to 12 
Gb for 8.4 million. Note that the RAM scales linearly with the number of unknowns (RAM 
oc MN) } whereas the CPU times scales with the 4/3 power of the number of unknowns (CPU 

4 

time oc It should be noted that the 8.4 million complex solution set consumed all of 

the RAM on the SGI ORIGIN 2000 and represents the largest number of complex unknowns 
that could be solved with this solver within the available RAM. 

In order to check the accuracy of the solver solutions, the authors use the reduction in 
the noise level from the entrance to the exit of the duct (i.e., AdB) as a metric. This is a 
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Number of unknowns, millions 
Fig. 3. Solver (VSS) statistics for a rigid-wall duct. 

physically more meaningful quantity than the error norm of the computed solution because 
it is the quantity that is perceived by the human ear as the noise source propagates down 
the duct. This metric has units of decibels, and is defined as 


AdB 


101og 10 


m o) 

E(L) 


(13) 


f H 1 

E(x) = -R e{p(x,y)u*(x,y)}dy (14) 

Jo Z 

where the superscript asterisk denotes the complex conjugate, Re{} denotes the real part of 
the bracketed quantity, and u is the normal component of acoustic particle velocity at axial 
location x. For the zero-flow calculations considered here, the acoustic velocity is related to 
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the pressure gradient via the axial momentum equation 

= (15) 

It can be shown from the exact solution that no sound is attenuated in a rigid-wall duct 
(i.e., AdB = 0). The solution vector {$} obtained from the equation solver is used to 
compute the noise level, AdB, numerically. This metric is then used to access the accuracy 
of the solver solutions. These statistics are plotted as a function of the number of complex 
unknowns in Figure 4. Note that noise level predictions computed from the equation solver 



Number of unknowns, millions 

Fig. 4. Noise level predictions from the solver solution vector. 

solution vector are in excellent agreement with the exact value of O.OdB. 

Statistics have also been computed for a soft-wall duct with identical dimensions as that 
of the rigid-wall duct, but with a wall impedance of ^(.r) = 0.5 — 0.5?’. Thq GPU time and 
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RAM requirements were minimally affected by the presence of the liner. The statistics 
for the soft-wall duct are therefore not presented, since they were nearly identical to those 
given in Figure 2. All computations presented here were computed by setting the loop level 
parameter [7] to unity. This value of the loop level parameter was observed to be optimal 
for these zero-flow computations. It has been observed that the choice of the loop level 
parameter could have a significant impact on the performance of the solver. For example, 
when the statistics given in Figures 2 and 3 were computed with a loop level of 6 (i.e. , the 
default loop level for the solver), both the CPU time and RAM requirements increased by a 
factor of 2. 

7 Conclusions 

The commercial version of the NASA General-Purpose Solver has been exercised to solve 
several problems in nacelle aeroacoustics that did not contain flow. Statistics show that the 
solver is capable of solving 10 thousand complex unknowns in as little as 0.01 sec and 8.4 
million complex unknowns in slightly less than 3 hours. The 8.4 million complex equation set 
represents an upper limit on the problem size that could be retained in memory (RAM) on 
an SGI ORIGIN 2000 with 12 Gb available RAM. The 8.4 million complex equation solution 
set also represents an upper limit of what would be expected in simulations of real laboratory 
experiments involving impedance eductions and is the largest aeroacoustics problem solved 
to date with this solution technique. This study was unable to detect any noticeable error 
in the solution, since noise levels predicted from the equation solver’s solution vector (i.e., 
with as many as 8.4 million complex unknowns) are in excellent agreement with the decibel 
levels computed from the exact solution. Statistics for rigid-wall ducts show that the solver 
RAM requirements and CPU times scale with the first and 4/3 power, respectively, of the 
number of unknowns. The performance of the solvers is minimally affected by the presence 
of the liner. 

Results presented in [3] indicate that the matrix equation described by Equation (7) must 


14 



be solved 20 to 40 times for successful impedance eductions at the low end of the impedance 
spectrum (i.e., frequencies below 3000 Hz). Although the solver allows the computations to 
be extended to the high frequency end of the impedance spectrum where 8 million equations 
need to be solved, impedance eduction still appears to be impractical on a single processor 
because a single pass through the solver requires nearly 3 hours of central processing unit 
time. This research, therefore, supports a recommendation that efforts be made to exploit 
the multi-processor capability of the solver so that aeroacoustic optimization studies become 
practical for large-scale aeroacoustics computations. This recommendation is in concert with 
the NASA 256-processor SGI ORIGIN 2000 system (to be upgraded to 512 processors) of 
the same type used in this study. 
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